Back to Article
Research Site Locations
Download Source

Research Site Locations

Alex Koiter

In [1]:
suppressPackageStartupMessages({library(tidyverse)
library(ggspatial)
library(sf)
library(rnaturalearth)
library(ggrepel)
library(raster)
library(patchwork)
})
Warning: package 'ggplot2' was built under R version 4.4.3Warning: package 'ggspatial' was built under R version 4.4.2Warning: package 'rnaturalearth' was built under R version 4.4.2Warning: package 'patchwork' was built under R version 4.4.2

Bounding box of region and sites

In [2]:
box = c(xmin = 455200, ymin = 5616000, xmax = 466200, ymax = 5626500)

ag_site <- data.frame(x = c(464124, 463318),
                         y = c(5622932, 5622121)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 3158) %>% 
  st_bbox() %>% 
  st_as_sfc()

forest_site <- data.frame(x = c(460822, 460189),
                         y = c(5619887, 5619206)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 3158) %>% 
  st_bbox() %>% 
  st_as_sfc()

locations <- data.frame(name = c("Forest", "Agriculture"), 
                        long = c(-99.560041, -99.513808), 
                        lat = c( 50.726076,  50.752731)) %>%
  st_as_sf(coords = c("long", "lat"),  crs = 4326)

  #st_transform(crs = st_crs(land))

cities <- data.frame(name = "McCreary", long = -99.493, lat = 50.775) %>%
  st_as_sf(coords = c("long", "lat"),  crs = 4326)

Land use data of region

Get data

In [3]:
minnedosa_data <- "lcv_minnedosa_2004_2006_shp"
if (file.exists(minnedosa_data)) {
  print("The data already exists!")
} else {
  download.file("https://mli.gov.mb.ca/landuse/shp_zip_files/lcv_minnedosa_2004_2006_shp.zip", 
                destfile = "lcv_minnedosa_2004_2006_shp.zip")
  unzip("lcv_minnedosa_2004_2006_shp.zip", exdir = "lcv_minnedosa_2004_2006_shp")
}
[1] "The data already exists!"

Summarize and crop

In [4]:
land <- st_read(here::here("./notebooks/lcv_minnedosa_2004_2006_shp/lcv_minnedosa_2004_2006.shp")) %>%
  st_crop(st_bbox(box)) %>%
  st_transform(4269) %>%
  mutate(Land_use = recode(CLASSNAME, "Agricultural Field" = "Agriculture", "Agri - Forage Field" = "Agriculture", "Water Body" = "Water", "Wetland - Treed Bog" = "Wetland", "Wetland - Marsh" = "Wetland", "Decidious Forest" = "Forest", "Open Decidious Forest" = "Forest", "Mixedwood Forest" = "Forest", "Forest Cut Blocks" = "Forest", "Coniferous Forest" = "Forest","Roads Trails Rail Lines" = "Other", "Sand and Gravel" = "Other", "Cultural Features" = "Other", "Range and Grassland" = "Forage")) %>%
  st_make_valid() %>%
  st_simplify(dTolerance = 0.0001, preserveTopology = F) %>%
  group_by(Land_use) %>%
  summarise(do_union = TRUE) %>%
  mutate(Land_use = factor(Land_use, levels=c("Agriculture", "Forest", "Forage", "Water", "Wetland", "Other")))
Reading layer `lcv_minnedosa_2004_2006' from data source 
  `C:\Users\koitera\Dropbox\Spatial variability\spatial-variability-soil-manuscript\notebooks\lcv_minnedosa_2004_2006_shp\lcv_minnedosa_2004_2006.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 160653 features and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 429972 ymin: 5542070 xmax: 549972 ymax: 5656880
Projected CRS: UTM
Warning: attribute variables are assumed to be spatially constant throughout
all geometriesWarning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
argument preserveTopology cannot be set to FALSE when working with ellipsoidal
coordinates since the algorithm behind st_simplify always preserves topological
relationships

DEM data

Freely available from https://maps.canada.ca/czs/index-en.html Data type is Elevation. You provide the clipping area and then submit the data request

In [5]:
topo <- raster::aggregate(raster(here::here("./notebooks/DEM.tif")), fact = 3) %>%
  crop(extent(land)) %>%
  as.data.frame(xy = T)

Hydro data

Freely available from https://maps.canada.ca/czs/index-en.html Data type is the CanVec. You provide the clipping area and then submit the data request

In [6]:
water_linear <- st_read(here::here("./notebooks/canvec_190514_89905/water_linear_flow_1.shp")) %>%
  st_transform(crs = st_crs(land)) %>%
  st_make_valid() %>%
  st_crop(st_bbox(land))
Reading layer `water_linear_flow_1' from data source 
  `C:\Users\koitera\Dropbox\Spatial variability\spatial-variability-soil-manuscript\notebooks\canvec_190514_89905\water_linear_flow_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 5160 features and 44 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -99.9082 ymin: 50.54151 xmax: -99.12695 ymax: 51.16603
Geodetic CRS:  GCS_North_American_1983_CSRS98
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Canada wide map

In [7]:
canada <- ne_states(country = "Canada", returnclass = "sf") %>%
  st_transform(crs = 3348)

locations2 <- data.frame(name = c("Center"), long = c(-99.989671), lat = c(50.742554)) %>%
  st_as_sf(coords = c("long", "lat"),  crs = 4326) %>%
  st_transform(crs = st_crs(land))

Maps

Canada

In [8]:
p1 <- ggplot() +
  theme_bw() +
  layer_spatial(canada, fill = "white") +
  layer_spatial(locations2, shape = 23, colour = "red", fill = "red", size = 2) +
  annotation_scale(location = "bl",
                   height = unit(0.05, "cm")) +
  annotation_north_arrow(location = "br", 
                         height = unit(0.5, "cm"),
                         width = unit(0.5, "cm")) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm"))

Land use

In [9]:
p2 <- ggplot() +
  theme_bw(base_size = 10) +
  geom_sf(data = land, aes(fill = Land_use)) +
  scale_fill_manual(values=c("#458B00", "#006400", "#CDAD00", "#009ACD", "#CDB38B",  "#C1CDCD"), name = "Land use") +
  geom_sf(data = water_linear, colour = "#009ACD", lwd = 0.5) +
  geom_label_repel(
    data = head(locations),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.01, 
    size = 3) +
  geom_label_repel(
    data = head(locations),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.01, 
    size = 3) +
  geom_sf(data = ag_site, fill = NA, colour = "red", lwd = 1) +
  geom_sf(data = forest_site, fill = NA, colour = "red", lwd = 1) +
  geom_label_repel(
    data = head(cities),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.01,
    size = 2) +
  annotation_scale(location = "bl",
                   height = unit(0.05, "cm")) +
  annotation_north_arrow(location = "br", 
                         height = unit(0.5, "cm"),
                         width = unit(0.5, "cm")) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.key.size = unit(5, 'mm'),
        legend.box.spacing = unit(0, "pt")) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0))

DEM

In [10]:
p3 <- ggplot() +
  theme_bw(base_size = 10) +
  geom_raster(data = topo , aes(x = x, y = y, fill = DEM)) + 
  scale_fill_viridis_c(name = "Elev. (m)  ", limits = c(282, 800), breaks = c(300, 550, 800)) +
  geom_label_repel(
    data = head(locations),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.01,
    size = 3) +
  geom_sf(data = ag_site, fill = NA, colour = "red", lwd = 1) +
  geom_sf(data = forest_site, fill = NA, colour = "red", lwd = 1) +
  geom_label_repel(
    data = head(cities),
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.01,
    size = 3) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  annotation_scale(location = "bl",
                   height = unit(0.05, "cm")) +
  annotation_north_arrow(location = "br", 
                         height = unit(0.5, "cm"),
                         width = unit(0.5, "cm")) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm"),
        legend.position = "bottom",
        legend.box.spacing = unit(0, "pt")) 
In [11]:
p6 <- p1 + guide_area() + p2 + p3 + plot_layout(ncol = 2, guides = "collect") + plot_annotation(tag_levels = 'a', tag_suffix = ')')
p6
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Figure 1: a) Map showing the location of the study sites within Canada. Location of the two study sites and nearby town of McCreary, and regional b) land use, and c) topography.